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Modern extensions of density functional theory such as the density functional theory plus U and 
the density functional theory plus dynamical mean-field theory require choices, including selection 
of variable (charge vs spin density) for the density functional and specification of the correlated 
subspace. This paper examines these issues in the context of the “plus U” extensions of density 
functional theory, in which additional correlations on specified correlated orbitals are treated using 
a Hartree-Fock approximation. Differences between using charge-only or spin-density-dependent 
exchange-correlation functionals and between Wannier and projector-based definitions of the corre¬ 
lated orbitals are considered on the formal level and in the context of the structural energetics of the 
rare earth nickelates. It is demonstrated that theories based on spin-dependent exchange-correlation 
functionals can lead to large and in some cases unphysical effective on-site exchange couplings. Wan¬ 
nier and projector-based definitions of the correlated orbitals lead to similar behavior near ambient 
pressure, but substantial differences are observed at large pressures. Implications for other beyond 
density functional methods such as the combination of density functional and dynamical mean field 
theory are discussed. 


I. INTRODUCTION 

Modern theories of electronic structure can be formally 
constructed in terms of functionals of observables of in¬ 
terest whose stationary points deliver the values of the 
observables— Practical use of this formal construction 
requires a choice of variables and of approximations to 
the functional. Perhaps the most common choice is den¬ 
sity functional theory^^ (DFT), which can be formulated 
as an effective action that is a functional only of the elec¬ 
tron density^ (ie. not spin resolved). While DFT is in 
principle exact, existing approximations have had diffi¬ 
culty capturing phenomena related to the formation of 
local magnetic moments. For example, neither the local 
density approximation^ (LDA) nor the generalized gra¬ 
dient approximation (GGA)£ provide correct accounts of 
the structural energetics of layered and spinel manganites 
that exhibit cooperative Jahn-Teller distortions associ¬ 
ated with the high spin state of Mn 3+ because at ambi¬ 
ent pressure both LDA and GGA incorrectly predict that 
the Mn ion is in a nominal \t\ g e° g ) low-spin configuration 
instead of the proper high-spin \t\ g e 3 ) state. Indeed it 
seems intuitively clear that it would be prohibitively dif¬ 
ficult to construct a functional based only on the density 
that could capture this sort of effect, while a functional 
of the spin-density might have a robust approximation 
which can capture this physics. Additionally, a func¬ 
tional of the spin-density will clearly allow predictions 
to be made about magnetism, and this avenue has been 
pursued since the inception of DFT— i 

Functionals of both the charge and spin density such 
as the the local spin-density approximation (LSDA)V— 
and the spin-dependent generalized gradient approxima¬ 
tion (SGGA)^il have been constructed. Such theories 
are often refered to as density functional theories, but in 


this paper we strictly distinguish terms, using the term 
density functional theory (DFT) to refer to theories such 
as the LDA and GGA that are based on a functional of 
the density only, and referring to theories such as the 
local spin density approximation (LSDA) or the spin- 
dependent GGA (SGGA) as spin-density functional the¬ 
ories (SDFT). 

SDFT theories perform far better than DFT theo¬ 
ries in describing the energetics of magnetic insulators, 
resolving, for example, the problems with manganites 
noted above— SDFT theories additionally make predic¬ 
tions about spin magnitudes and the nature of ordered 
states. However, the known implementations of SDFT 
fail to correctly describe many aspects of the physics and 
structure of strongly correlated electron systems, for ex¬ 
ample providing qualitatively incorrect structures for the 
rare-earth nickelate o 12 ’ 13 (see, e.g. Ref. HI for additional 
examples). 

These difficulties motivated the construction of new ef¬ 
fective action theories that depend not only on the den¬ 
sity or the spin-density, but also on additional proper¬ 
ties of a subspace of orbitals for which correlations are 
believed to be relevant-^. Subspaces which have been 
treated in this way include the transition metal d orbitals 
in transition metal oxides and the lanthanide/actinide 
/ levels in heavy fermion compounds. Various differ¬ 
ent variables can be defined from the subspace of corre¬ 
lated orbitals. In this paper we focus on the historically 
important and currently widely used choice of the site- 
local spin and orbitally resolved density matrix associ¬ 
ated with the correlated subspacei^; however, we expect 
that our findings are relevant to other variable choices; 
in particular to the case of dynamical mean field the¬ 
ory where the additional variables are the components of 
the local Green’s function. A straightforward functional 
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to use with the local spin-resolved density matrix is the 
Hartree-Fock energy functional, defined using site-local 
matrix elements of the Coulomb interaction. The func¬ 
tional resulting from this combination of a local Hartree- 
Fock functional of a correlated subspace and a standard 
functional of the density or spin-density is commonly re¬ 
ferred to as a “+U” extension of density functional the¬ 
ory. 

The paper that introduced the “+U” approach^ em¬ 
ployed a functional of the electron density only, so is 
referred to here as a DFT+U approach. However, the 
vast majority of subsequent papers^ employ SDFT func¬ 
tionals and are referred to here as SDFT+U approaches. 
Despite the many successes of the +U methodology, ba¬ 
sic points including the rationale for choosing SDFT+U 
in preference to DFT+U and the factors influencing the 
construction of the correlated subspace have not been 
clearly discussed. While the general formalism may be 
applied with any choice of orbital, it is important in prac¬ 
tice to choose orbitals that are optimal in the context of 
the other approximations used in constructing the theory. 
While there is no clear prescription for doing this, the lo¬ 
cal nature of the approximations which will ultimately 
be used suggests that it is sensible to choose the corre¬ 
lated subspace to consist of well localized ‘atomic-like’ 
orbitals. These orbitals may be constructed from by pro¬ 
jecting onto a set of localized orbitals defined to lie within 
atomic spheres. This choice is natural, given that many 
basis sets for electronic structure already utilize projec¬ 
tors. Projectors are used in various beyond-DFT meth¬ 
ods including (S)DFT+Ui& as well as DFT+DMFT— 
Alternatively, Wannier functions may be used to con¬ 
struct the correlated orbital sets for beyond DFT cal¬ 
culations. Various forms of the Wannier function have 
been used for DFT+U2^ and DFT+DMFT— includ¬ 
ing the projected Wannier function, iVth-order muffin-tin 
orbitals^, and the maximally localized Wannier function 
(MLWF)^-35. 

In this paper we describe the physical differences be¬ 
tween DFT+U and SDFT+U and provide guidelines to 
enable researchers to choose between them. We further 
compare the effect of different correlated orbitals sets 
(Projector vs Wannier) on energy calculations within the 
(S)DFT+U method. We also present a comprehensive 
discussion of the issues arising when the +U methodolo¬ 
gies are combined with the Projector Augmented Plane 
Wave (PAW) formalism widely used to perform efficient 
(S)DFT calculations. In addition to the formalism we 
provide a quantitative application in the context of the 
relation between crystal structure and energetics of the 
rare-earth nickelates. This family of materials provides 
a useful benchmark because its members exhibit a struc¬ 
tural phase transition which is not correctly captured 
either by DFT or by SDFT calculations. We com¬ 
pare DFT+U and SDFT+U, using both projectors and 
Wannier functions to construct the correlated subspace. 
These results can also be directly compared to our recent 
DFT+DMFT total energy calculations for the same class 


of materials— 

We note in passing that an additional important issue 
in DFT+U and SDFT+U theories is the so-called double 
counting correction, introduced to account for the fact 
that the local interactions denoted by U and J are to 
some degree present already in the (S)DFT. The double 
counting issue has been addressed in great detail in previ¬ 
ous wort— and is not critical to the issues examined here. 
Therefore in this paper we use the conventional definition 
of the “fully localized-limit” double-countin g 36 ’ 37 . 

The rest of this paper is organized as follows. Section 
El presents the basic formalism. Section Ell provides 
a careful discussion of the issues involved in combining 
the +U formalism with the projector augmented plane 
wave method. Section HVI nresents expressions for forces, 
needed in optimizing structures. Section [V] compares 
the DFT+U and SDFT+U (with projector and Wan¬ 
nier definitions of the correlated subspace) predictions for 
the structural properties of the rare earth nickelates as a 
function of unit cell volume, while section IVIIII provides 
a comparison of predictions for the rare earth nickelate 
phase diagram and the equilibrium volume. Section IIXI 
contains a summary and conclusions. Relevant compu¬ 
tational details are given in Appendix [X] 

II. FORMALISM 

In this section we present explicit formulae for the en¬ 
ergy functional, using a variant of presentation of the 
DFT+DMFT functional given in Ref. El We derive the 
total energy functional for SDFT+U, and then obtain 
the DFT+U functional as a special case. 

The SDFT+U total energy functional is defined by 

E[p°,n™] = Tr[(H^)] - Tr[V£ xc • p°] 

- Tr[V£ t ■ n Ta ] + E Hxc [p a } + E mt [n TrT } (1) 

where p a denotes the charge density of electrons with spin 
a (we neglect spin-orbit coupling here for simplicity) and 
n Ta is a density matrix within the correlated subspace of 
an atom r. Here the bracket () means that the eigenstates 
of Hy are summed over for the eigenvalues less than the 
Fermi energy. 

The functional E Hxc is the familiar Hartree and 
exchange-correlation energy functional of the DFT the¬ 
ory to be used and V§ xc = 5E Hxc /8p G is the correspond¬ 
ing Hartree-exchange-correlation potential. 

E mt is the combination of a Hartree-Fock potential 
energy E pot defined within the correlated subspace and 
a double counting correction E DC introduced to remove 
from this potential the parts already included in the un¬ 
derlying DFT: 

E int = E pot _ e dc ( 2 ) 

In the applications presented here we follow the com¬ 
mon practice in the literature by choosing the correlated 
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subspace to be particular orbitals m of spin a on partic¬ 
ular atoms r in the solid; for example the 3 d orbitals in 
a first row transition metal ion or the 4/ orbitals in a 
lanthanide ion, so that 

rppot _ \ \rrcri<J2criCT2 Tjcricr2&2cn X )n Tr71 n T ‘ 72 

m i m 3 rn^m^ 7721771377147712 ' u o’io '2 )' L m\m 2 

T,m 1 m 2 

m 3 m 4 °T <7 2 

( 3 ) 

where the are the site-local matrix elements 

of the Coulomb interaction within the correlated sub¬ 
space, appropriately renormalized by the solid-state en¬ 
vironment. This paper will present an application of the 
formalism to the case of transition metal d-orbitals where 
(if spin orbit coupling is neglected), the may 

be parametrized by the Slater integrals F°, F 2 , and 
F 4 . We typically further assume that the d-wave func¬ 
tions are sufficiently similar to their free-space forms that 
F 4 = 0.625 F 2 . It is then conventional to define on-site 
interaction U and the Hund’s coupling J via F° = U, 
F 2 = (14/1.625) 

The double counting energy E DC is the subject of con¬ 
tinuing discussion in the literature but these subtleties 
are not relevant here. We adopt the ‘fully localized limit’ 
form given by 

E DC m = jN t (N t - 1 ) - ^ N r( N r - 1 )> ( 4 ) 

a 

where N° is the total number of electrons of spin a on site 
r and N t = is the total on-site charge density in 

the correlated subspace. Note that use of SDFT (in other 
words the dependence of E Hxc \p a ] on the spin density) 
means that that the double counting potential E^ c in 
Eq.[4]also depends on spin indices. 

The corresponding interaction potential V^ t is ob¬ 
tained as V? t = dE int /dn Ta : 

Knt = \ T ’ m ’ a )( V pot,mm' ~ V5c)( r > m '> ( 5 ) 

m,m' 


nuclei at the position R (the interaction between nuclei is 
not explicitly denoted but is included in the formalism). 

The physical state of the system at zero temperature is 
obtained by extremizing Eq.[l]with respect to variations 
in p a and n Ta . To perform the extremization, we solve 
the eigenvalue problem of Hy and find the eigenstates 
|'I , f k ) (i is a band index and k is a momentum in the 
first Brillouin zone). We then determine the local charge 

density p a { r) = Yik fik l 4 'ifc( r )| 2 (/ is the Fermi function 
and it is evaluated at zero temperature (T=0)). The on¬ 
site density matrix is also determined from |\H) using for 
example the method described in Section IIII Bl or Section 
1III Cl Finally, we require consistency between the p a and 
n Ta and the Kohn-Sham and V m t potentials they imply. 

Once self-consistent solutions of p a and n Ta are ob¬ 
tained, the total ground-state energy E can be obtained 
from the value of Eq.[l]at the stationary point of both p a 
and n Ta . It is also useful to cast the total energy func¬ 
tional in Eq(T]is a slightly different form, both for analysis 
and for technical reasons specially when using Wannier 
functions. The SDFT+U total energy functional can be 
decomposed into the SDFT energy ( E SDFT ), the KS en¬ 
ergy correction ( E AKS ), and the interaction energy cor¬ 
rection ( E mt ) (defined only within the correlated sub¬ 
space) as follows: 

E = E sdft + E aks + E int , (10) 

where 

E sdft = Tr [(H^s)} - Tr[VZ xc ■ p°] + E H ™[p°] (11) 
and 

e aks = Tr [( R-)] - Tr[Vf nt ■ n - TV[<i% s >]. (12) 

DFT+U is a special case of the SDFT+U in which the 
exchange-correlation energy depends only on the total 
density p, i.e., $Ks(p a ) -+ ^ks(p)- However, spin de¬ 
pendence is retained in the correlated subspace, so n rcr 
is still considered to be spin dependent. Thus the total 
energy functional in Eq.[T] becomes 


where 

F/ot.mm' = ^ ] ( Um 1 mm 2 m ' — U mirnm ' m2 • 6cr 1 a)n m ^ m2 

( 6 ) 

and 

VZh = U(N d -±)-J(N2-±). ( 7 ) 

Finally, the effective Hamiltonian H^j is 

Hu = Hks + Fj nt > (8) 

where the Kohn-Sham Hamiltonian H^g is 

H a KS [p\ R] = ~V 2 + Vext [R] + VSJjT]. (9) 

H°j is the analogue of the Kohn-Sham Hamiltonian H^g 
of SDFT. V ex t is an external potential arising from atomic 


E[ Pl n T ”] = Tr[(H?j)} - Tr[V Hxc • p] 

- Tr[VP t ■ n Ta } + E Hxc [p\ + E int [n Ta ], (13) 

where 

HZ = H K s + V? t . (14) 

The rest of the formalism carries through as before, ex¬ 
cept that the exchange-correlation potential now depends 
only on p and therefore the double counting correction is 
taken to be spin-independent: 

E dc [N t \ = ^N t (N t -1)- ^N T (N T - 2), (15) 

implying 

V DC [N T \ = U{N T - l -)-i{N T -l). (16) 

Thus in DFT+U theories spin dependence arises only 
from the properties of the correlated subspace (which af¬ 
fect the rest of the system via hybridization). 
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III. CORRELATED ORBITALS AND THE 

IMPLEMENTION OF (S)DFT+U USING THE 
PROJECTOR-AUGMENTED PLANE WAVE 
METHOD 

A. Overview 

In this paper, the DFT portion of the functional 
lEci. 1111) is implemented using the projector augmented 
wave (PAW) method^. The main idea of PAW is to cir¬ 
cumvent treating the computationally inconvenient core 
states by use of a linear transformation which relates an 
all-electron wavefunction |i/>) to a pseudo wavefunction 
\ip). The transformation requires the addition of augmen¬ 
tation terms which can be expanded using a projector 
function | p) and the resulting KS Hamiltonian contains 
additional terms arising from the augmentations, but 
because the resulting pseudo-wavefunction is smoothly 
varying computations are much more efficient. 

In the PAW method, the SDFT energy functional can 
be split into three terms: 

e sdft _ ^ ~ Zc ] + E 1 [n 1 , n Zc ] - E 1 [n 1 , n, n Zc ] 

( 17 ) 

where E is the pseudo-energy term, E 1 is the on-site all¬ 
electron energy term, and E 1 is the on-site pseudo-energy 
term, h is the pseudo charge density, n 1 is the on-site all¬ 
electron charge density, n 1 is the on-site pseudo charge 
density, n is the compensation charge between n 1 and h 1 
such that n 1 + n has the exact same moment as n 1 , h Zc 
is the pseudized core density, and n Zc is the all-electron 
core density. 

The effective s for generating the pseudo wavefunc- 
tion \ip) is now given by extremizing E SDFT in Eci.fTTl 
with respect to the pseudized charge p a : 

Hks = ^ 2 +Veff+J2 ( 18 ) 

where u e // is the effective pseudo one-particle potential 
obtained using v e ff = dE a ^ —• Dij, D \-, and Djj are po¬ 
tentials conjugate to the density matrix of the augmenta¬ 
tion part Pij , ie, D ZJ = Dh = |g, and Dh = |g. 

Eq.[T7]can be cast into a similar form as Eq.fTTl 

E SDFT^ p a ] = ^f ik $ ik \H£ s \ip ik )+E% c AW {h,n,n 1 ,n 1 }. 

ik 

( 19 ) 

The PAW double counting correction E F ^ W also con¬ 
tains a pseudo part and an augmentation part: 

E d c AU [n, h, n 1 , fi 1 ] = E dc [n, n\ + E l dc [n l } - E^ c [n 1 ,n]. 

( 20 ) 

The PAW double counting correction E F ^ W should not 
be confused with the double-counting correction E DC re¬ 
quired in the interaction functional. The derivation of 
the above equations and the explanation of each term 
are given in Ref. HI. 


The treatment of the +U interactions in the PAW for¬ 
malism depends on the prescription used to construct the 
correlated subspace. Accordingly, the remainder of this 
section is divided into two parts, one dealing with the 
projector formalism ('subsection IIII Bl) and one with the 
Wannier formalism (subsection IIII Cl) . 

B. Projectors: ortho-normalization 

We begin with the projector method, in which the com¬ 
ponents n™ n of the correlated orbital density matrix n 
appearing in Eq.[l]are obtained by projecting the Kohn- 
Sham wavefunction ip onto the spherical harmonics Yj m 
inside an atomic sphere centered on atom r with the ra¬ 
dius r£, i.e., 

= (21) 

ik 

Here i is a band index and k is a wavevector in the first 
Brillouin zone, /ik is the Fermi function evaluated at 
T=0 throughout our paper and P^ nn is the projector func¬ 
tion on atom r defined by 

(r'\PUr) = - <)© (r T < r T c ) (22) 

where r T = r — R T is the position vector defined with 
respect to the atomic center R T and 0(:r) is the step 
function such that 0(x) = 1 if x < 0 and 0(x) = 0 if 
x > 0.— Note that if the ferrni function is removed from 
Eq.EUrmd the sum is taken over all bands i and momenta 
k then standard completeness relations imply that 

Omn = I Ann I Vw) ~ O n 5 mn (23) 

ik 

is a diagonal matrix, whose normalization depends on the 
choice rl of sphere cutoff. 

Within the PAW formalism, n™ n is computed from the 
pseudo-wavefunctions \ip) and a pseudo-projector P Ta as 

<,n = £/*$«kfcfe)- (24) 

ik 

The pseudo-projector is defined in terms of an appropri¬ 
ate set | (p a ) of solutions to the Schroedinger equation for 
a reference atom in free space as 

P™n = 'E\P“)( < t , a\PU<t>b){p b \ (25) 

ab 

where the |p) are the PAW projector functions conjugate 
to the d> n . In practice, we use the implementation in the 
VASP code. 

The wave functions defined by the PAW projector pro¬ 
cess do not constitute an orthonormal set because the 
sum is only over a subset of states and a choice of sphere 
radius is made. Thus the overlap matrix 

0™ = £<^J^> (26) 

ik 
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is neither diagonal nor possessing correctly normalized 
eigenvalue. To obtain a properly ortho-normalized den¬ 
sity matrix n within the correlated subspace we define 

TCn = E (° TCT ) w • <?n' • (onj! 2 - (27) 

m'n' 

This procedure of the ortho-normalization of the projec¬ 
tor function is used in some DFT+DMFT implementa¬ 
tions^. 

Within the PAW formalism, the effective Hamilto¬ 
nian Hy can be obtained by varying the energy func¬ 
tional Eq.[T]with respect to the pseudized charge p a = 

X) ik /ik I V’ik ) (V'ik I : 


TTCT _ TT(7 

H U — h ks 


dE d 


dp° 


= H 


dn T 


KS 


dp a 


(28) 


The dr ' d ~ term in Eq. [28] is difficult to evaluate because 
the overlap matrix O in Eq. [27] also varies implicitly due 
to the change of |V'ik) - For simplicity, we assume that 
the ortho-normalization effect is fully incorporated in the 
electronic V^ t [n Ta ] term while the change of the density 
matrix via p is computed from the un-normalized n Ta : 

dfi Ta dn Ta 

. — <29, 

ij 

This interaction potential part is expanded with the basis 
of the projector |p) and it can be added to the augmenta¬ 
tion part of H^s * n Eq.EDl Therefore, the Hamiltonian 
Hfj is given by 

Hu = ~ 2^ 2 Veff + E I Pi)(Dij + Djj — b\j + Vij)(pj\ 

(30) 

where 


V l3 =V° nt [n™]-{&\P T \<t>i)- ( 31 ) 


The total energy functional within PAW can be ob¬ 
tained as follows: 

E[p a 1 n Ta ] =^f ik (^ ik \H^ik) + E^ c AW [n,h,n 1 ,n 1 } 

ik 

-Tr(V° t [n T<r }-n T n + E M [n™]. (32) 

The interaction energy correction E mt (E pot (Eq. [3]) - 
E dc (Eq. ]4[) ) term and the interaction potential correc¬ 
tion V^ t term are computed using the orthonormalized 
density matrix n^ n . In practice, the band energy correc¬ 
tion term Tr (V° lt [n Ta ] ■ n Trj ) is computed only updating 
the V^ t term while the density n Ta is obtained from un¬ 
normalized projector functions. In this way properly or¬ 
thonormalized correlated orbitals can be used with only 
a slight modification of the PAW formalism. 


C. MLWF orbitals 

Here, we derive the +U formalism in the case where 
the correlated subspace n™ n is defined by Wannier func¬ 
tions. We follow the approach used in our previous anal¬ 
ysis of the DFT+DMFT formalism—. In this subsection 
we present the formalism purely for DFT+U, as all of 
the comparisons between projectors and Wannier in this 
study will take place in the context of DFT+U. The gen¬ 
eralization to SDFT+U is however straightforward. 

Wannier functions are discussed at length in the 
liter ature^r— . Here we make only a few remarks. First, 
the construction of a Wannier functions requires the spec¬ 
ification of a hybridization window W, a range of ener¬ 
gies from which the states used in the construction of the 
Wannier functions are defined. This energy range should 
encompass both the correlated orbitals and the orbitals 
which directly hybridize with the correlated ones. For ex¬ 
ample, in the case of the rare-earth nickelates, Wannier 
functions are constructed from an energy window (« 11 
eV wide) including the full Ni-3d and 0-2 p manifolds. 
By construction the Wannier functions provide a com¬ 
plete orthonormal basis for states within the energy win¬ 
dow so it is not necessary to introduce an overlap matrix. 
A continuous infinity of choices of Wannier basis exists; 
here we choose a “maximally localized” (MLWF—) basis 
set that minimizes the sum of Wannier function spreads 
((r 2 )-(r) 2 ) and also perform an additional orbital rota¬ 
tion as described below Eq. [33] We denote the resulting 
states as |W Rt ) where r labels an atom within a unit 
cell, R denotes a lattice vector, and n is an orbital index. 

The projection of the DFT Hamiltonian onto the 
MLWF basis set is 

H^n T ' = {Wn' T '\H KS \W^)- (33) 

As discussed e.g. in Ref. [li|, for the Wannier functions 
pertaining to the correlated states we perform a rotation 
in the orbital indices to minimize the off-diagonal terms 
of the on-site correlated-state }j^- corr ’ T ’ Ilcorr - T j n the mn 
subspace. 

The DFT+U calculation solves the eigenvalue prob¬ 
lem of the Hks + matrix within the hybridization 
window W. One should note that the V^ t term is spin- 
dependent while Hks has no explicit spin dependence 
and has parameters determined by the total density. 

The density matrix rj within the hybridization win¬ 
dow, which includes the correlated subspace as a subset, 
is obtained from the eigenvalues and eigenfunctions of 
Hks + y^nt as 

V™ n = E (34) 

l£W. k 

The density matrix n™ n in the correlated subspace is 
a sub-block of r}™ n . Our basis choice in the nm space 
means that the off-diagonal terms are negligible; n™ n ~ 

$nm ■ 
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The band energy correction E AKS is then given by 

E AKS = Tr(H KS -r 1 )-Tr(H KS -r ] 0 ), (35) 

where the rf^ n is computed via Eq. [M] but using 
the eigenfunctions and eigenvalues of Hks rather than 
Hks + V° nt , 

The interaction energy can be defined using Eq.[3]with 
the calculated density matrix 1 Eq J 34I) . In our application 
to the nickelates we construct this term using the Slater- 
Kanamori Hamiltonian as defined in our previous paper— 
and for ease of reference present the results in the same 
notation. In Ref. [l3| the on-site intra-orbital interaction 
is given as u, the Hund’s coupling is j the inter-orbital 
interaction is u — 2 j and the exchange and pair-hopping 
terms do not contribute in the Hartree-Fock approxima¬ 
tion used here. The parameters are related to the U and 
J defined elsewhere in the paper by u = U + (8/7) J and 
j = (5/7)J. For the specific application to the nicke¬ 
lates we took u=6.14eV and /=0.71eV (corresponding to 
U= 5eV and J=leV). 

The interaction potential energy is then 


formalism, therefore one can adopt the the same calcula¬ 
tion already implemented in VASP. In the case of Wan- 
nier functions, additional implementation is needed to 
compute the derivative of Wannier function \W^) with 
respect to R. 

The force functional following from the projector cor¬ 
related orbital set can be derived from the energy func¬ 
tional in Eg . 1321 taking derivatives with respect to R pro¬ 
duces the same functional form as the PAW force func¬ 
tional (see e.g. Ref. for details) except that the eigen¬ 
values e n and eigenfunctions ip n are obtained by solv¬ 
ing Hy. The force has terms due to the change of the 
pseudized core charge density fiz c via the explicit move¬ 
ment of the ionic positions, the change of the compen¬ 
sation charge fi itself, and the change of the projector 
functions | p) as the ions are moved: 


pDFT+U 



d(Hu - e n (l + J2ij \Pi)Qij(Pj \)) 
<9R 



F 1 

dhzc 

+ F 2 

' dfi' 

+ F 3 

r<9|p}<pli 


<9R 


Lc»rJ 


<9R 


(38) 


E pot = u 1 


, T W^4 . 


0-2 j) 

m^m' ,r 


t t rJ, 
n rn Urn' 


+ (« - 3 j) 


(36) 


m>m / ,rcj 


while the double counting energy is 


where 


F 3 


' d\p)ip\ - 
. <9R 


- 520a + D]j - D]j + Vij - e n qij ) 

n,ij 


x fni^r 


, 9\Pi)(Pj 
' <9R 


!*")■ 


(39) 


E DC = ^N d (N d - 1) - ^N d {N d - 2), (37) 

where here Na is the total occupancy of the d levels on a 
Ni ion. 


IV. THE FORCE FUNCTIONAL 

The force functional is defined in terms of the deriva¬ 
tives of the energy functional with respect to atomic po¬ 
sitions. Here we present the force functional correspond¬ 
ing to the DFT+U version of the energy functional in 
Eq. [T3J the SDFT+U forces can be derived similarly. 
The specifics depend on the formalism. As before, the 
PAW formalism as implemented in VASP is utilized and 
we present forces for the projector, while we outline the 
differences for the case of a Wannier based correlated 
subspace. 

The DFT+U force functional consists of two parts, 
namely the same functional form used in DFT except 
that the DFT Fermi function is replaced by the DFT+U 
density matrix utilizing DFT+U eigenvalues and eigen¬ 
functions and an additional force term derived from the 
E mt energy term. The computation of the forces requires 
consideration of the derivatives of the correlated orbital 
density matrix with respect to R, which in turn arise 
from the changes of the correlated orbital wave functions 
as the atomic positions change. The derivative of the 
projector orbital is already computed within DFT force 


qij is the correction to the overlap matrix given by 
((f>i\(f>j) — The explicit expressions of F 1 and 

F 2 are the same as DFT forces given in Ref. The 
implicit changes of the Hamiltonian Hfj via the density 
n, n 1 , n 1 , and fi are always cancelled out exactly against 
the change of terms in Ea.l32l 

The evaluation of a force term from E mt requires the 
derivative of the correlated orbital density matrix, i.e. 

and the result depends on the choice of correlated 
orbital sets. Within the projector scheme, the force term 
arising from the implicit change of interaction energy cor¬ 
rection E mt fEa. [32l : via the ortho-normalized density 
matrix n Trj is given by 


dE int \n Ta } dn Ta 
dn Ta dR 



(40) 


Here, the calculation of term is complicated by the 
R-dependence of the overlap matrix O (Ea ETl) . In practi¬ 
cal applications the derivative of the density matrix with 
respect to the ionic position R is thus approximated to 
the change of the un-normalized n T<J via the derivative of 
the projector function which is already present in PAW 
force (Eq. [39]): 


dn Ta 

F^czVUn™}'^. ( 41 ) 

Taking a derivative of the — Tr(V mt [h Ta ] ■ n Ta ) term in 
Eq. [35] with respect to R leads to a term — U mt [h T£T ] • 
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ri ” R which cancels out the term in Eq. [41] and a term 
— dv " 1 • n r<T which cancels out the implicit change of 


V mt [n Tl 7 ] term in H^. 

The force functional using Wannier functions can be 
derived in a similar way as the projector functions except 
that the term needs to be computed explicitly in 

terms of the change of Wannier functions with R. How¬ 
ever, we have not yet implemented this. Instead, we have 
determined the minimum energy structure in the phase 
space of pressure and bond length difference as defined 
in our previous paper—. 


V. APPLICATION TO STRUCTURAL 
PROPERTIES OF RARE-EARTH NICKELATES 

A. Overview 

In following sections we investigate the general issues of 
interest in this paper, namely different correlated orbital 
sets (Projector vs Wannier) and background electronic 
structure methods (DFT+U vs SDFT+U), in the spe¬ 
cific context of the structural properties of the rare-earth 
nickelates, I 7 N 1 O 3 . In these materials the basic struc¬ 
tural motif is the NiOg octahedron. At some values of 
temperature, pressure and R, the materials exhibit a uni¬ 
form phase in which all octahedra have approximately 
the same mean Ni-0 bond length. In other parameter 
regimes the materials exhibit a two-sublattice dispropor- 
tionated phase in which the octahedra on one sublattice 
have a mean Ni-0 bond length ~ 0.1 A shorter than the 
octahedra on the other sublattice. The materials are im¬ 
portant for the present study because this basic struc¬ 
tural property is closely linked to a fundamental elec¬ 
tronic property, namely whether the material is a metal 
or a correlation-driven site-selective Mott insulato r 32 ! 41 . 
This linkage means that obtaining a correct description 
of the structural properties poses a critical test for the 
electronic structure methods. 

Our previous studies^^ showed that (S)DFT+U does 
not provide a quantitatively accurate description of the 
experimental structural and metal-insulator transition 
phase diagram of I?Ni 03 series. A particular difficulty is 
the prediction that the ambient-pressure ground state of 
LaNiOs is bond-disproportionated and insulating when 
the actual material has a non-disproportionated i?3c 
structure and is metallic. A closely related deficiency of 
the DFT+U method is an overestimation of the critical 
pressure of the metal-insulator transition for materials 
where the ambient pressure ground state is insulating. 
DFT+DMFT methods produced much better results. 
It is also the case that DFT+U (and DFT+DMFT) 
wrongly predict that the ground state is ferromagnetic. 

However the trends found in the DFT+U calculations 
were found to track the trends found in the DFT+DMFT 
calculations. For example the DFT+U T —> 0 structural 
phase boundary in the pressure-tolerance factor plane 


was offset by a certain pressume from the DFT+DMFT 
phase boundary, indicating that the difficulty is simply 
that the DFT+U methods overestimate (to a consider¬ 
able degree) the stability of the insulating state. Thus 
since the aim of the current investigation is understand 
how different formulations of the theory affect basic is¬ 
sues of structure and energetics, rather than to accurately 
model material properties, we can use the +U approxi¬ 
mation as a flexible and inexpensive computational labo¬ 
ratory with the expectation that the similarity of trends 
noted above suggests that the general findings will be 
applicable to DFT+DMFT as well. 

We present results computed as described in Ap¬ 
pendix El In the current paper, we use the Vienna 
Ab-initio Simulation Package (VASP) 39 ! 42 which adopts 
the PAW formalism. For the exchange-correlation DFT 
functional, we use a generalized gradient approxima¬ 
tion (GGA) with the Perdue-Burke-Ernzerhof (PBE) 
functional^ and also adopt a local density approxima¬ 
tion (LDA) if necessary. We take the correlated orbitals 
to be atomic-like Ni-centered d orbitals using the projec¬ 
tor method and assume the additional interactions have 
the form given in Eq. [3] Unless otherwise specified, we 
use U= 5eV and J= leV for all computations. 

We present results for three members of the material 
family: LuNiOs (strong insulator at ambient pressure), 
NdNiOs (insulating but near the phase boundary at am¬ 
bient pressure), and LaNiOs (metallic at ambient pres¬ 
sure) . In the rest of this section we introduce the materi¬ 
als and compare the DFT+U and SDFT+U results using 
both the GGA and LDA functionals for bond lengths and 
energetics. In the following section we discuss the pro¬ 
jector vs Wannier issue and in a third section present 
implications for computed phase diagrams. 


B. Bond disproportionation vs volume 

We used the VASP implementation of (S)DFT+U us¬ 
ing the GGA functional to perform full structural relax¬ 
ations of the three compounds from the low-symmetry 
bond-disproportionated structure with the unit cell vol¬ 
ume constrained to take particular values (see Fig.[l]). 
The correlated orbital basis set was treated using the 
same orthonormalized projector for all calculations (note 
that the orthonormalization of the basis set required that 
we modify the VASP energy and force formulas accord¬ 
ing to Eq. EH- In many cases a disproportionated struc¬ 
ture with two inequivalent NiOg octahedra was found; 
for these cases we computed the difference 5a in mean 
Ni-0 bond length between Ni sites on different sublat¬ 
tices. Generically if 5a ^ 0 the band structure exhibits a 
gap at the fermi level (for very small disproportionation 
amplitudes the gapping may not be complete), but for 
simplicity of presentation we do not consider the elec¬ 
tronic structure here. 

Results are shown in Fig. |1] At positive compression 
(V — V 0 < 0) the differences between SDFT+U and 





GGA Funtional 




Figure 1. (Color online) Average Ni-O bond-length differ¬ 
ence 5a as a function of volume computed using the GGA 
functional with DFT+U (top panel) and SDFT+U (lower 
panel) for LuNiOa (red square), NdNiCU (blue diamond), 
and LaNiOa (green circle). U=5eV and J=leV are used. 
The same orthonormalized projector method is used to con¬ 
struct the correlated subspace in all cases. Vo is the zero- 
pressure volume computed for the given compound by the 
given method. 


DFT+U are quantitative, but large. We see that con¬ 
sistently across the material family SDFT+U predicts 
a higher 5a at given compression and similarly predicts 
that a higher critical compression is needed to drive the 
structural transition (5a —> 0) than does DFT+U. How¬ 
ever, at negative compression (V — Vo > 0) the difference 
is qualitative: DFT+U predicts a monotonic increase of 
5a values as (V — Vo)/Vo increases while SDFT+U cal¬ 
culations indicate a reduction of 5a as the cell volume is 
increased and ultimately a reentrant structural transition 
(seen in the data for LaNiOa and expected for the other 
materials from the downward curvature). 

Fig. [2] shows the effect of varying the Hunds coupling 
J on the computed bond disproportionation for LaNiOa. 
In the DFT+U calculations, increasing J from 0 to 1 eV 
has a dramatic effect, while a further increase to 2 eV has 
a weaker effect, suggesting a saturation as J is increased. 
On the other hand, the SDFT+U results show almost 
no J-dependence, indicating that the spin dependent ex¬ 
change potential in SDFT already effectively includes a 
large on-site J and suggesting that J is not needed when 
performing SDFT+U calculations. This could be prob¬ 
lematic for SDFT+DMFT calculations, where the dy¬ 
namical effect of J are typically important. 

The upper portions of the two panels of Fig. [2]further 
show that within DFT+U the choice of DFT method 
(LDA vs GGA) produces quantitative but not quali¬ 
tative differences over the volume range investigated, 
with in particular the LDA+U exhibiting a smaller 5a 
at given J and volume, consistent with the known ten- 



(V-V a )/V 0 [%] 


Figure 2. (Color online) Bond disproportionation 5a plot¬ 
ted against relative volume change (V — Vo )/Vo computed for 
LaNiOs for different Hunds coupling J at U= 5eV using both 
(a) GGA and (b) LDA functionals and different methods in¬ 
dicated in the legends (DFT+U vs SDFT+U). The same or¬ 
thonormalized projector method is used to construct the cor¬ 
related subspace in all cases. Vo is the zero-pressure volume 
computed using the given method and J=leV. The change 
of Vo at different J is negligible for SDFT+U and rather no¬ 
table for DFT+U (see Fig.[4]) but it does not affect any results 
discussed here. 


dency of the PBE GGA functional used here to over¬ 
estimate magnetism^. Alternatively, in SDFT+U the 
choice of DFT method produces a qualitative difference, 
with the SGGA+U method indicating reentrance of the 
non-disproportionated phase at a small positive relative 
volume while no indication of reentrance is found in 
the LSDA+U calculations. Interestingly, The J= 2eV 
GGA+U calculations also suggest that reentrance of the 
undistorted phase would occur at larger relative volumes, 
consistent with the notion that the SDFT methods imply 
a large (perhaps excessively large) J already at the SDFT 
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level. Taken together these results also suggest that the 
mathematical origin of the reentrant transition is a (pre¬ 
sumably unphysical) effect of large J. The disproportion- 
ated phase may be understood as a hybridization density 
wave corresponding to relatively strong Ni-0 bonding at 
the short-(5 — a site — so although the precise connection 
is not clear at this point we may speculate that the reen¬ 
trance is related to unphysically large spin-dependence of 
level shifts of the Ni-d relative to O -p states, weakening 
the Ni-0 singlet bond that produces the distortion. 


netic moment on J, with magnetic moment increasing 
with J with the dependence becoming weaker as the sat¬ 
uration value M = 1 /is is reached and the critical vol¬ 
ume for the magnetic transition also being J-dependent. 
The difference between LDA+U and GGA +U results 
reflects the stronger tendency toward magnetism charac¬ 
teristic of the PBE-GGA functional. In effect, PBE-GGA 
already contains a certain degree of local exchange. In 
contrast, SDFT+U calculations in both GGA and LDA 
produce large moments at all volumes, and with negligi¬ 
ble J dependence. 


GGA Funtional 




(V-V 0 )/V 0 [%] 


Figure 3. (Color online) The magnetic moments per Ni atom 
in LaNiC >3 obtained using different J as a function of pressure 
using both (a) GGA and (b) LDA functionals and different 
methods indicated in the legends (DFT+U vs SDFT+U). The 
same orthonormalized projector method is used to construct 
the correlated subspace in all cases. Vo is the zero-pressure 
volume computed using the given method and J=leV. 


The interplay between the DFT functional (LDA vs 
GGA), the value of Hund’s J and the physics of the dis¬ 
proportionation instability are also evident in the study 
of the magnetic moments presented in Fig. [3j The 
DFT+U results reveal the expected dependence of mag¬ 


VI. ENERGETICS 



(V-v 0 )/v 0 [%] 


Figure 4. (Color online) The total energy of LaNiOs as a func¬ 
tion of relative volume difference computed using GGA+U 
(top) and SGGA+U (bottom) with an orthonormalized pro¬ 
jector definition of the correlated orbitals at J values indi¬ 
cated. At each volume and for each method the structure is 
relaxed and the energy of the relaxed structure is presented. 
The zero of energy is chosen at a compression of 5% in all 
curves. 


In this section, we compare the DFT and SDFT predic¬ 
tions for energies. We restrict attention to the GGA and 
SGGA density functionals, define the correlated states 
via orthonormalized projectors, and focus on LaNiOs. 
For each relative volume, the structure was relaxed and 
then the energy was evaluated. Fig [3] displays the de¬ 
pendence of the total energy on the normalized volume 
difference for different J values. 

As found in the previous section’s analysis of the dis¬ 
proportionation amplitude and magnetic moment, sub¬ 
stantial differences between DFT+U and SDFT+U are 
found. The DFT+U energy curve depends substantially 
on J, changing rapidly as J is increased from zero and 
saturating as J becomes large. Remarkably even the 
equilibrium volume is J-dependent. The SDFT+U en¬ 
ergy has negligible J dependence and is similar to the 
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J= 2eV DFT+U result again suggesting that the SDFT 
exchange correlation functional in effect contains a J 
which (for the PBE-GGA case studied here) is substan¬ 
tially larger than the J ~ leV values believed to be phys¬ 
ically reasonable. 
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Figure 5. (Color online) Contributions to DFT+U energy 
functional (cf Eq. 1101 ) computed for LaNiCU as a function of 
relative volume, a) The DFT contribution and b) the corre¬ 
lation correction CC contribution ( E AKS + E mt ) are decom¬ 
posed from Fig. [4] using the same relaxed structure at each 
volume and compared for different methods (DFT+U (the 
top panel) vs SDFT+U (the bottom panel)) and different J 
values. 


Fig. [5] presents a decomposition of the energy into 
the DFT contribution (E DFT , panel a) and the corre¬ 
lation correction (CC) contribution (E AKs + E mt , lower 
panel), as defined in Eq. [lOl for DFT+U (upper half 
of each panel) and SDFT+U (lower half of each panel). 
The DFT term E DFT contains the structural contribu¬ 
tion while the CC term expresses the correlation physics. 
E dft is not monotonic in unit cell volume, expressing 
the basic physics of chemical bonding. Alternatively, the 
CC contribution decreases monotonically as the volume 
is increased for both DFT+U and SDFT+U and for all J 


values, expressing the enhancement of correlation occur¬ 
ring when hybridization is decreased and showing that 
the equilibrium volume predicted by the correlated cal¬ 
culation is larger than that predicted from the DFT con¬ 
tribution alone. 

Fig. [5] indicates that for the SDFT+U method neither 
E dft nor the CC term has significant J dependence 
because the SDFT method already includes a large lo¬ 
cal exchange contribution. Alternatively, in the DFT+U 
calculation both terms have some J dependence. The 
DFT+U CC term changes dramatically as J is increased 
from 0 to leU, accounting for the noticeable change of 
DFT+U energetics from J=0 to leV in Fig. [4] but does 
not change much as J is further increased because the 
moment is saturated (see Fig. [3]). It is also interesting to 
note that the DFT+U CC energy at J > leV is com¬ 
parable to the SDFT+U CC energy at J = 0, further 
confirming the large value of J implicit in the SDFT 
method. In the DFT+U method, some dependence of 
E dft on J occurs as J is increased from leU to 2eV, 
and it is interesting to note that this is the J range where 
suggestions of reentrance are visible in the DFT+U cal¬ 
culation. This behavior again indicates that the unusual 
reentrant behavior is related to a rearrangement of the 
band-structure by an unphysically large Hunds coupling. 



Figure 6. (Color online) The occupancy of the correlated 
orbitals expressed as the number of d electrons Nd per Ni 
atom computed for LaNiCU at J values indicated for GGA+U 
(the top panel) and SGGA+U (the bottom panel). 

As extensively discussed elsewhe re 45 ’ 46 , the occupancy 
Nd of the correlated orbitals provides useful insights 
into the physics of strongly interacting electron systems. 
Fig. [HI displays the d occupancies computed for LaNiCU 
for the parameters whose energies are shown in Fig. [5] 
Within DFT+U, increasing the Hunds coupling J leads 
to a decrease in Nd , a signal of stronger correlation arising 
from effectively smaller p-d hybridization. In SDFT+U 
the correlations in this sense are already stronger at the 
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SDFT level (ie. SDFT+U yields smaller N ( j). and adding 
additional J does not change the situation. 


VII. CHOICE OF CORRELATED ORBITAL 



(V-V 0 )/V 0 [%] 


Figure 7. (Color online) The average Ni-0 bond length dif¬ 
ference 5a as a function of reduced volume computed for ma¬ 
terials indicated using DFT+I7 as implemented with the en¬ 
ergy functional of Eq.[lO] Both maximally localized Wannier 
functions (filled symbols, solid lines) and ortho-normalized 
projectors (open symbols, dashed lines)) are compared. In¬ 
teraction parameters of U= 5eV and J=leV are used for 
the projector-based calcualtions while the equivalent values 
u=6.14eV and j=0.71eV (Slater-Kanamori parameterization) 
are correspondingly used for the Wannier construction. 

In this section we study the effect of the choice of cor¬ 
related orbital on the calculated results. The DFT+U 
method is used; SDFT+U is not considered in this sec¬ 
tion. Fig[7] presents the bond disproportionation ampli¬ 
tude 5a versus reduced volume for different rare earth 
nickelates using either MLWF (filled symbols, solid lines) 
or ortho-normalized projectors (open symbols, dashed 
lines) for the correlated subspace. The qualitative trends 
of 5a as a function of reduced volume are similar for 
both correlated orbitals. Substantial differences appear 
only for very large compression, where the Wannier ap¬ 
proach enhances the tendency to the bond disproportion- 
ated states, though only for NdNiC >3 and LaNiC> 3 . The 
origin of the difference is not clear at present, but may 
have to do with the fact that size of the Wannier func¬ 
tion varies as the volume and the bond length difference 
change, while the projectors are defined using a fixed ra¬ 
dius. 

We now turn to the effect of orthonormalization, com¬ 
paring in Fig. [5] structural relaxation calculations per¬ 
formed using the orthonormalized projector orbital to 
calculations and performed using the unnormalized pro¬ 
jector implemented in VASP. Normalization has a par¬ 
ticularly important effect in the small volume region of 



(V-V 0 )/V 0 [%] 


Figure 8. (Color online) Average Ni-0 bond length difference 
5a graph as a function of volumes obtained using DFT+U. 
The normalization effect of the correlated orbital is inves¬ 
tigated by comparing both the ortho-normalized projector 
correlated orbital (filled symbols, solid lines) and the un¬ 
normalized projector (open symbols, dashed lines). LuNiC >3 
(red square), NdNiC >3 (green diamond), and LaNiC >3 (blue 
circle) results are displayed for comparison. 

L 11 MO 3 , where the critical pressure calculated using the 
unnormalized projector is overestimated and probably in¬ 
correct. However, normalization has no noticeable effect 
on the structural relaxation of LaNiCU. As the volume is 
expanded the consequences of normalization are seen to 
be minor. This discrepancy in the small volume region 
may arise from an overestimate of the spectral weight in¬ 
side the Ni atomic sphere, leading to a mis-estimate of 
the correlation energy. 


VIII. THE PHASE DIAGRAM OF ANIO3 AND 
THE EQUILIBRIUM VOLUME V 0 

In this section, we compute the structural phase di¬ 
grams and equilibrium volume Vo of the rare-earth nick¬ 
elates 7 ?NiC >3 obtained using DFT+U and SDFT+U as 
functions of the reduced volume for different R ions. 
MLWF and project definitions of the correlated orbitals 
are also explored. 

Fig-El displays the structural transition phase diagram 
of rare-earth nickelates .RMO 3 in the plane of reduced 
volume (V — Vo)/Vo) and choice R of rare earth ion 
computed using the different beyond DFT methods dis¬ 
cussed in this paper. The structural transition is be¬ 
tween the P2i/n structure {5a >0) and Pbnm struc¬ 
ture (da=0) for all i?Ni 03 except rhombohedral LaNiC> 3 . 
For LaNiCU, the transition associated with the bond- 
disproportionation separates the R3 structure {5a >0) 
and the R3c structure (6a=0). 

All DFT+U and SDFT+U results produce critical 
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Structural Phase Diagram 

- DFT+U (MLWF) 


0 SDFT+U (norm-Projector) 
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Figure 9. (Color online) Structural phase diagram of the rare- 
earth nickelates as a function of the volume compression and 
rare earth ion, computed using DFT+U (square dots) and 
SDFT+U (circle dots). The un-normalized projector (open 
dots) and ortho-normalized projector (filled dots) are com¬ 
pared for both methods. The MLWF correlated orbital result 
(diamond dot) is also shown for DFT+U, yielding some sim¬ 
ilarity to the phase boundary to the ortho-normalized pro¬ 
jector. The experimental data (black dash-dot lines) are also 
given for comparison^ (see also Ref. 03 and 0 ). 


compression for the transition which is too large relative 
to experiment. (DFT+DMFT produces results in much 
better agreement with experimen t 12 1 13 ). The SDFT+U 
method (red circular dots) exhibits more rapid variation 
of the critical compression with change of R ions than 
does DFT+U (purple square dots). 

The effect of the ortho-normalization in a correlated 
orbital varies depending on the functional. SDFT+U im¬ 
plemented using the un-normalized projector as adopted 
in VASPfcircular open dots and dashed lines, also shown 
in Ref. Ill) moderately reduces the critical pressures (fa¬ 
voring the Sa >0 region) compared to the same SDFT+U 
implemented using the ortho-normalized projector (cir¬ 
cular filled dots and solid lines). In contrast, DFT+U 
with the ortho-normalized projector (square filled dots 
and solid lines) more substantially moves the critical line 
toward the structural phase with da=0 ( Pbnm ) for the 
heavy rare earths. 

DFT+U implemented using the MLWF basis set (blue 
diamond dot and solid line, also shown in Ref. [H with 
u=5eV and j'=leV of Slater-Kanamori parametrization) 
also produces a similar phase boundary to the orthonor- 
malized projector DFT+U calculation compared to other 
methods. Therefore, we deduce that Wannier and or- 
thonormalized projectors yield similar behavior. 

The compression (vertical axis) of the structural phase 
diagram in Fig. [9] is defined as the relative change of vol¬ 
ume compared to the equilibrium volume Vo computed 
within each theoretical method. Vq is determined as the 



Figure 10. (Color online) Equilibrium volume Vo calculated 
using SGGA (filled symbols) and LSDA (open symbols) for 
-RNiCU with R= La, Sm, Nd and Lu using DFT+U (square 
dots) and SDFT+U. For SDFT+U, ortho-normalized projec¬ 
tor (green diamond) and un-normalized projector (red circle) 
calculations are compared. The experimental volumes at the 
ambient pressure are depicted by a black dashed line with 
open pentagonal dots. LDA(+U) results are shown only for 
LaNiOs since the pseudo potentials for other rare earth ions 
except La are not available. 

volume of the minimum energy from an energy vs volume 
curve and the atomic positions at each volume are ob¬ 
tained by minimizing the inter-atomic forces. Results are 
displayed in Fig. [TUI First, we discuss results obtained us¬ 
ing pure SDFT within the LSDA and SGGA approxima¬ 
tions. The calculated Vo with GGA exchange-correlation 
functional (filled symbols) is larger than the experimen¬ 
tal one (dashed line), while the Vo obtained with the 
LDA exchange-correlation functional (open symbols) is 
smaller than the experimental value. This behavior is 
well known in the DFT literature. 

We next observe that the calculated Vo values com¬ 
puted using SDFT+U are rather sensitive to the ortho¬ 
normalization of correlated orbitals. The orthonormal- 
ized projectors lead to substantially larger equilibrium 
volumes than the un-normalized projectors, leading to 
better agreement with experiment for LDA based func¬ 
tionals and worse agreement for GGA, though in neither 
method is the agreement particularly good. 


IX. CONCLUSION 

We studied different formulations of DFT+U and 
SDFT+U in the context of total energy calculations in 
the rare-earth nickelates. The correlated subspace was 
constructed in three different ways: maximally local¬ 
ized Wannier functions, orthonormalized projectors, and 
non-orthonormalized projectors. We computed the Ni-0 
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bond length difference Sa as a function of pressure, the 
structural phase diagram describing the transition be¬ 
tween the bond-disproportionated structure (Sa >0) and 
the no bond-disproportionated structure (Sa= 0) as func¬ 
tions of pressure and the rare-earth ions, and also the 
equilibrium volume. 

SDFT+U and DFT+U show qualitatively different 
behavior in some circumstances. In particular, the 
SGGA+U results of the structural transition in the rare- 
earth nickelates show a re-entrant transition with pres¬ 
sure, and this is not observed in GGA+U calculations 
that are performed with a reasonable on-site exchange 
J = leV. However, increasing J to 2eV in GGA+U 
produced similar results to SGGA+U, implying that the 
SGGA spin-polarized exchange correlation functional re¬ 
sults in a large effective on-site exchange. SDFT+U 
based on the LSDA exchange-correlation functional re¬ 
sults in a more reasonable effective J, meaning that 
LSDA+U results using <7=0 are similar to LDA+U with 
J ~leV. The reentrant transition at negative pressure 
does not occur within DFT+U calculations using the 
LDA functional for J <2eV. Our results suggest that 
there is no need to use an on-site exchange when per¬ 
forming SDFT+U (ie. set J = 0), and this is effectively 
equivalent to using the approach of Dudarev et al^. Ad¬ 
ditionally, our results imply that the SGGA should be 
used with caution given its overemphasis of local ex¬ 
change. More generally, DFT+U with an appropriately 
chosen J can largely recover the qualitative behavior of 
SDFT+U. 

We demonstrated that orthonormalized projectors be¬ 
haved rather similarly to MLWF near ambient pres¬ 
sure; although notable differences are evident for NdNiCU 
and LaNiOa under large pressures. Additionally, the 
un-normalized projector, as implemented in VASP, can 
lead to notably different results, especially at high pres¬ 
sures. Within SDFT+U, the equilibrium volumes are 
substantially increased when computed using the ortho¬ 
normalized orbitals compared to un-normalized orbitals. 

Given that (S)DFT+U is equivalent to 
(S)DFT+DMFT when the DMFT quantum impu¬ 
rity problem is solved within Hartree-Fock, we expect 
the general findings of this study to be applicable 
to (S)DFT+DMFT as well. Given our finding that 
DFT+U with an appropriately chosen J can largely 
recover the qualitative behavior of SDFT+U, our work 


supports the long held tradition of basing dynamical 
mean field extensions on DFT theories rather than 
SDFT theories. This is particularly important in cases 
where the dynamical effects of J are crucial. Our 
results have broad implications for the application of 
SDFT+U, given that in the nickelates the choice of 
SDFT functional leads to dramatic differences in the 
effective on-site exhcange interaction. A very similar 
conclusion was reached in a study of the spin crossover 
molecule Fe(phen) 2 (NCS)z^. 
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Appendix A: Computational details 

The (S)DFT+U formalism is implemented in VASP 
using the projector functions | p) as a correlated orbital 
set. The Hamiltonian for (S)DFT+U is given by Ea.l30l 
and the total energy and force equations are also derived 
in Eg .[T+jl and EaJ551 However, the VASP implementation 
adopt the un-normalized density matrix n Ta as Ea.l24l In 
the current paper, we compare the VASP implementation 
to the ortho-normalization of n Ta to give rise to n Ta as 
derived in E ciITTl Fortunately, the VASP implementation 
provides both DFT+U (LDAUTYPE=4) and SDFT+U 
(LDAUTYPE=1) methods. 

For performing the summation of k points in the Bril- 
louin zone, we used the tetrahedron method^. When 
using projector correlated orbitals, a k-point mesh of 
6x6x6 (for the Pbnm and P2i/n structures) and 
8x8x8 (for the LaNiOa i?3c and R3 structures) are 
used with an energy cutoff of 600eV. When using the 
Wannier functions as correlated orbitals, k points meshes 
of size 10 x 10 x 10 is used for Pbnm and P2i/n, while 
16 x 16 x 16 for R3c and R3. 
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